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I. INTRODUCTION 


Guided missiles can be classified into four categories 
depending on launch and target position characteristics. The 
categories are Air-to-Air, Air-to-Surface, Surface-to-Air, and 
Surface-to-Surface missiles. 

Another classification among missiles is the guidance 
system of the missile. The missile can be command or homing 
guidance. 

In the command guidance system the missile and target are 
continuously tracked and guided from one or more friendly 
vantage points, and the necessary path for intercept is 
computed and relayed to the missile. 

In the homing guidance system, the missile has a homing 
device onboard which can detect the target and gives the 
necessary path directions for intercept to the missile. The 
homing missile is further subdivided into classes having 
active, semiactive, and passive guidance systems. Active 
detection is when the missile illuminates the target, i.e., 
with a radar, and receives the reflected signals. Semiactive 
detection is when the target is illuminated by a source other 
than the missile and the missile receives the reflected 
Signals. Passive detection is used when the target is the 
source of energy, and the missile detects signals that 


propagate from the target. 


Each of the missiles in the above categories will employ 
one or more of the three guidance laws. These laws are Pursuit 
Guidance, Line-of-Sight Guidance, and Proportional Guidance. 
The first portion of the missile flight path may use one of 
the guidance laws but the terminal phase of flight may be best 
suited for another. 

The present work addresses the design and evaluation of a 
semiactive Surface-to-Air missile using Proportional 
Navigation as the guidance law. A ground based target tracker 
will also developed with the target deviations in position and 
velocity relayed to the missile. 

Chapter II presents a description and comparison of the 
three different guidance laws. The Proportional Navigation 
guidance law will also be developed. In Chapter III the 
missile and target flight path models will be developed using 
the concepts of Chapter II and computer simulation studies 
will be performed. Chapter IV consists of a Luenberger 
observer design, and an evaluation of the estimator, and the 
guidance law over a range of conditions will be conducted. 
Chapter V consists of the development of the ground target 
tracker using the theory of the Kalman Filter and again, 
computer simulations will be included to determine the 
accuracy of the target tracker. 

All computer simulations are developed and conducted using 


the Matrix Laboratory (MATLAB) language. 


II. MISSILE GUIDANCE 


A. GUIDANCE LAW SELECTION 

The selection of a guidance law is a pre-requisite for 
determining the initial calculations for the model. The 
missile guidance system measures the error between the 
missile’s actual and desired course, computes the corrections 
necessary to reduce the error based on the guidance law 
selected, and gives commands to the autopilot to activate the 
controls required to achieve acceptable intercept of the 
target. The miss distance and the acceleration required by the 
missile are functions of the guidance law. 

1. Pursuit Guidance 

The pursuit guidance law is illustrated in Figure 1 

and is described as having the missile velocity vector 
directed toward the target at all times. The missile is always 
heading along the line-of-sight from the missile to the 
target. This guidance law is effective against slow moving 
targets, but the missile may lack sufficient maneuverability 
against fast moving maneuverable targets. 

2. Line-of-8ight Guidance 

| Line-of-sight guidance is used in a beam-rider type 
missile and is illustrated in Figure 2. This guidance law 


requires that the missile remain on a line joining the target 


TARGET 


9, * angle of missile heading 


B * ongte of line of sight 


MISSILE HORIZONTAL 





Figure 1. Pursuit Guidance Trajectory 


and the target tracker. The purpose of the target tracker is 
to maintain the antenna boresight pointing at the center of 
the reflecting area of the target. This guidance scheme 
normally requires a dedicated fire control system from launch 
to intercept (Ref.1]. 
3. Proportional Navigation Guidance 

This guidance law requires that the missile travel in 

such a way that its own rate of turn is proportional to the 
rate of turn of the line-of-sight from the missile to the 
IR, Figure 3 illustrates the proportional guidance scheme 
in which the rate of change of the missile heading is made 
proportional to the rate of change of the line-of-sight 


between the missile and the target. The fixed or variable 
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Figure 2. Line-of-Sight Trajectory 


multiple between the missile rate of turn and the rate of turn 
of the line-of-sight is called the navigation ratio (NR). The 
proportional navigation guidance law attempts to generate an 
acceleration command perpendicular to the line-of-sight. One 
way to achieve this could be lateral acceleration coupled with 
angular or angular rate commands to place the acceleration 
perpendicular to the line-of-sight. The advantage of this 
guidance law is in its effectiveness against maneuvering 
targets. Since proportional navigation guidance anticipates 
a kane te future position, it can attain a higher degree of 
responsiveness over the other two guidance laws. 


Figure 4 illustrates the proportional navigation 


scheme (Ref.2]. If the seeker head of the missile follows the 


TARGET 


l 
l 
l 
l 


I 


MISSILE HORIZONTAL 





Figure 3. Proportional Guidance Trajectory 


Ap = Ag = -£ (w RI + w R (25999) 


target, the target acceleration perpendicular to the line-of- 


sight will equal the acceleration of the R vector, where 


R - missile-target line-of-sight vector 

R = closing rate along R 

w = angular rate of change of R 

Ar = target acceleration perpendicular to R 
Ag = tangential acceleration of vector R 


The term (oR) represents the vectorial acceleration of 
R and the term d/do(oR) represents the rate of change of the 
tangential velocity. A missile acceleration, Ay, equal to the 


target acceleration A, at this point will make the line-of- 


TARGET 


MISSILE TARGET 
LINE OP SIGHT 





Figure 4. Proportional Navigation Scheme 


sight parallel to its original direction. Since the velocity 
R' is along the vector R, a missile/target intercept is 
assured. Therefore, from Equation (2.1), the executed missile 


acceleration commands should be 


Ay = 2 (WR) +(@R) (2.2) 


Since the direction of the velocity vector cannot be 
directly controlled, proportional navigation is achieved by 


‘controlling the commanded missile acceleration (acom). 


acom - V, Y (2.3) 


where 


Vu = the missile velocity 


4 the rate of change of the velocity vector 
Implementation of proportional navigation results in 


the following guidance law 


acom » NR V, Ó (2.4) 


where 
o = the rate of change of the line-of-sight 
With this definition, the equation for the rate of 


change of the velocity vector can be written as 


Y = NO (2.5) 


B. PROPORTIONAL NAVIGATION KINEMATICS 

From Figure 5 the following equations of motion are 
obtained. The three dimensional linear model will be developed 
in Chapter III, but for simplicity and ease of understanding, 
the fundamental equations will be first developed in the x and 
y planes. The velocity vector Vy is at an angle y from the 
established reference line. From the geometry of the problem, 


the missile flight path angle can be easily determined. 





Figure 5. Intercept Geometry for Proportional Navigation 


V, 
Y, = arctan ES (2.6) 
MX 


where Vix and Vy, are the components of the missile velocity 
vector in the x and y directions, respectively. The target 
flight path angle is found from the same geometry and is 


expressed as 


TY 


V. 


6 
í 


where Vu and Vr are the target velocities in the x and y 
directions, respectively. 


The missile-target line-of-sight vector is R. 


Rad MS (2.8) 


where X, and Y, are the x and y coordinates of the target and 
Xu and Y, are the x and y coordinates of the missile. 
The magnitude of the velocity vector for the missile and 


target (VQ) can be expressed as 


Vu = CC Vege 12 + (Va )? p? (2.9) 


VM UA (2.10) 


The line-of-sight angle (c) is defined as 


Y.- Y 
g- arctan | -—— —* (2 SE) 
Xr “Ay 


and the rotation rate of the line-of-sight (0) can be 


expressed as 


- Vrsintyp-0) -Vysin(Yy-0) (OR 
R 
where V, sin(y - 0) and Vy Sin(yy 7 c) are the target and 


missile velocity components normal to the line-of-sight. 


10 


III. KINEMATIC REPRESENTATION OF A SKID-TO-TURN MISSILE 


A. INTRODUCTION 

Two basic methods of controlling the attitude of a missile 
to achieve the acceleration commanded by the guidance law are 
skid-to-turn and bank-to-turn. In the skid-to-turn method the 
roll angle is held to a small quantity and is usually 
considered zero in initial calculations. The magnitude and 
orientation of the body acceleration vector is achieved by 
permitting the missile to develop both an angle of attack and 
a sideslip angle. The presence of the sideslip imparts a 
skidding motion to the missile. The bank-to-turn missile 
generally will develop higher lift accelerations than the 
skid-to-turn method, so that the missile requires banking 
maneuvers to properly direct the control vector. To achieve 
the desired orientation, the missile is rolled or banked so 
that the plane of maximum normal force is oriented in the 
desired direction. The present work assumes a skid-to-turn 


missile that is roll stabilized. 


B. ASSUMPTIONS 
A simplified, point mass representation of the missile 
equations of motion will be developed under the following 


assumptions. 
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The missile thrust exactly cancels drag. 


. The orientation of the missile can be described by the 


Euler angles that represent the flight path angles in 
pitch and yaw. 

The seeker head angle rate is a good estimate of the 
line-of-sight rate. 

For small angles of attack and sideslip, the velocity 
vector is assumed to be aligned with the body center 
line and if the missile speed is maintained nearly 
constant, the missile acceleration in the x direction 


is fairly close to zero. 


C. SIMPLIFIED EQUATIONS OF MOTION 


Figure 6 illustrates the flight path geometry of the 


missile or target, with the velocity vector aligned with the 


body centerline. From the flight path geometry, the missile 


and target velocity magnitudes are expressed as 


for the 


for the 
The 
defined 


ratios. 


Vu = E (Og)? (CV) * (V)? ] 7? 


(3.1) 
missile, and 
Ve tal Vee a IS (3.2) 
target. 
angles of attack (y pitch) and sideslip (y yaw) are 


as kinematic quantities depending only on velocity 


From Figure 6, the missile and target flight path 
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Figure 6. Flight Path Geometry 


angles can be defined as 


Yw pitch = arctan — o | (3.3) 
VT (Vaa) IP (Vu) )) 
Yx yaw = arctan ( a 
M Va (3.4) 
V. 
Vo BESE tray | on 
ET NÉS po 


13 


t "ry 
Yr yaw = arctan us (3.6) 


If the velocity vector is not aligned with the body 


centerline, the angle of attack is defined as 


V 
Yy pitch = arctan Ca 


7 (3.7) 


with the sideslip angle in the yaw plane remaining unchanged. 


Figure 7 illustrates the sightline geometry. 


target 
(Xr, Yr, Z4) 


missile 


(Xy, Yu, 24) 





Figure 7. Sightline Geometry 


From the sightline geometry, the line-of-sight angles can be 


determined. 
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o pitch = arctan 





zacz, | 


VI (Xp=Xu) 2#(Yp-Y.)2) EN 
Yr=Yy 
o yaw = arctan | —— 
Xr=Xy (3.9) 
The range vector (R) can be defined as 
R = A O IAS AA 5) (3.10) 


The velocity components can be found from the Euler angle 
transformation. The derivation of the Euler angles can be 
found in Reference 4. The x, y, and z directions are taken as 
the longitudinal, lateral, and normal axis of the missile, 
respectively. The corresponding angles that represent the 
angular displacements are $ (roll), 6 (pitch), and y (yaw). 
With the missile assumed to be roll stabilized, the roll angle 
will be taken as zero and not included. The Euler 
transformation matrix is given as 

x cos(0)cos(y) cos(8)sin(w) -sin(8)] |Z 
y| = -sin(y) cos (y) 0 AN 
z sin(0)cos(y) sin(0) sin(y) cos(8) 
where I points north, J east, and K down. The transformation 
matrix represents the total transformation from the inertial 


coordinate system to the missile body axes. 
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From the above matrix, the linear components of the 
missiles velocity can be found using the fight path angles in 


the yaw and pitch planes. 


(3.12) 
Vay = Vy(cos (8) sin(y) ) (3-13) 
Viz = Vy(-sin(0)) (3.14) 


A skid-to-turn missile is controlled by generating the 
required acceleration commands in the pitch and yaw plane. 
From Reference 2, the commanded acceleration in the pitch and 


yaw planes are defined as 


ACOM icen = Vu P1€£CH Yy pitch (3.15) 


Acom,,, = V, yaw Yy yaw (3.16) 


The velocity components of the missile in the pitch and yaw 


planes are defined as 


V, pitch » V, cos( Y, yaw - a yaw) (3508 


Vy Yaw 7 V, COS( Y, pitch) (3.18) 


The autopilot of the missile will convert the commanded 


accelerations into angular rates 7, pitch and 7, yaw. These 
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angular rates are passed to the control surface servos and are 
converted to the necessary fin deflections required to steer 


the missile to the desired course. 


D. SYSTEM SIGNAL FLOW GRAPH 

Figure 8 shows the total system signal flow graph in the 
y direction, where at launch the missile/target line-of-sight 
is taken as the x axis. The outputs from the seeker is the 
seeker head angle rate, which when multiplied by the 
navigation ratio, becomes the input to the autopilot. The 
autopilot provides the signals required for the missile 


dynamics. 


E. LATERAL AUTOPILOT 

As stated previously, the autopilot receives commands from 
a guidance computer and processes these commands to control 
the deflections of the control surface. There are three 
principle requirements when designing an autopilot. These are, 
quick response, stability, and robustness. The autopilot 
should be able to handle broad variations of the aerodynamic 
parameters. For example, the missile may be required to 
operate over a broad range of Mach numbers and at various 
altitudes which will effect the normal force coefficient, 
which is a function of the Mach number and the altitude. 

In this simulation, a simple lateral autopilot will be 


modeled as a first order lag with a time constant (T) of 0.33. 
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Figure 8. System Flow Graph 
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This has been demonstrated in past simulations to adequately 
approximate a more detailed autopilot. 

The lateral autopilot controls the motion of the missile 
in the pitch and yaw planes. As in the case of this study, a 
symmetrical cruciform missile, the pitch and yaw autopilots 
are often identical so only one will need to be modeled. With 
the assumption that the angle of attack is very small, the 
velocity of the missile is aligned with the missile 


centerline. From Figure 8, we have the following expression 
Yy= “KYy+ KNRB (3.19) 


where 8 is the seeker head gimble angle rate and 7, is the 
angular acceleration of the missile flight path. K is equal to 


( 1/T). The transfer function of the autopilot is given as 


Yaa (S) . K NR (3.20) 
B (S) S+K | 


F. SEEKER HEAD DEVELOPMENT 

A homing head, when mounted in an airborne missile, is 
called a seeker. The purpose of the seeker is to detect, 
acquire and track a target by sensing radiation or reflected 
energy from the target. In the present study, a narrow field- 
of-view seeker gimballed to the airframe will be designed. The 
seeker maintains the target within this narrow field-of-view 


by rotating the platforn. 
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If the platform is inertially stabilized, rotation is achieved 
by applying torques which are proportional to the target's 
displacement. From the centerline of the field-of-view, the 


equation of motion is 


T=IB (3.219 


where T is the applied torque, I is the moment of inertia, and 
8 is the angular acceleration of the seeker head angle. From 
Figure 8 and equation (3.21) the resulting equation of motion 
is 


B= + = -KGB -K,B + K, o (3:295 


where K, and K, are functions of the time constants of the 
seeker. The transfer function of the seeker head can then be 


expressed as 


B (s) _ K, 


o (S) S+KS+K, (3.23) 


G. CONTINUOUS TIME STATE EQUATIONS 
1. Missile Dynamics 


Given the continuous time state equation 


X (t) =Ax (t) + Bu (t) (3.24) 


the missile or target state equation is defined. 
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2. Autopilot State Equation 
The pitch and yaw autopilots are identical, so the 


following state equation for the autopilot is defined as 


23 T 
0 -3 


3. Seeker Head State Equation 


Ypi tch 


Yu Yaw 


(3.26) 








. ; l B pitch 








The seeker head continuous state equation with a time 


constant of 0.1 seconds in the pitch plane is defined as 


Be 





0 1 | B pitch 
-100 -20||B pitch 





0 
Sa o (3.27) 





The yaw autopilot is identical to the pitch autopilot, so the 
above equation can be easily determined in the yaw plane. 
4. Continuous to Discrete 
Given the continuous time state equations, the 


discrete time equations are defined. 
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x (k+ 1) ® x(k) + VT u(k) 


(3.28) 


y (k+ 1) =Cx(k) + Duík) 


(3.29) 


The simulation study will use the discrete time representation 


of the given continuous time equations. 


H. SIMULATION 

This section presents the results of the missile/target 
engagement using the proportional navigation scheme. The 
following assumptions are made: 

1. The missile is limited to 20 g's in either the yaw or 
pitch plane. 

2. The target is capable of a 5 g maneuver. 

3. The seeker head system is noise free. This will forma 
basis from which the following chapters can be referred 
to. 

4. The minimum difference between the targets position and 
the missiles position, will be used as an estimate of 
the miss distance parameter. 

5. The origin is taken as coordinates (0,0,0) in the x,y,z 
plane. 

1. Constant Velocity Target 
The first scenario will be a constant velocity target 


with the following initial conditions. 
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Vux(0) = 2000 feet/sec 
X, (0) = 20,000 feet 
Z (0) = 50 feet 


Vix(0) = 1000 feet/sec 
All other initial conditions are zero. As shown in Figure 9, 
a successful intercept occurred for the constant velocity 
target. The minimum distance between the target and missile 
was 0.7665 feet. Figure 10 is a plot of the missile/target 
trajectories looking down on the X-Y plane. It shows the 
target starting at 20,000 feet on the X axis with the missile 
starting at the origin. Figure 11 is a plot of the 
missile/target trajectories in the Y-Z plane with the missile 
starting at the origin. Figure 12 is a plot of the missile 
commanded acceleration parameter. The acceleration increases 
rapidly during the initial phase of the engagement as force is 
applied to the missile and then drops off until just prior to 
intercept. Figure 13 shows the seeker head angle rate in the 
yaw plane and Figure 14 is a plot of the seeker head angle 
rate in the pitch plane. These parameters also increase 
rapidly during the first few seconds as the target is starting 
out on its flight profile and then decreases in time as the 


missile is on course due to the acceleration commands. 
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Figure 9. Range vs Time 


MISSILE-TARGET TRAJECTORY IN THE X-Y PLANE 





Figure 10. Trajectory in the X-Y Plane 
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MISSILE- TARGET TRAJECTORY IN THE Y-Z PLANE 





Figure ll. Trajectory in the Y-Z Plane 





Figure 12. Commanded Acceleration 
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Figure 13. Seeker Head Angle Rate in the Yaw Plane 





Figure 14. 


Seeker Head Angle Rate in the Pitch Plane 
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2. Constant Acceleration Target 
The second scenario will be a constant acceleration 

target with the following initial conditions: 

Vux(0) = 2000 feet/sec 

Xr (0) = 20,000 feet 

Zr (0) = 50 feet 

Viy (0) = 500 feet/sec 

Ap (0) = -5.0*32.2 feet/sec^2 

Ay (0) = 3.5*32.2 feet/sec”2 

Ar (0) = 3.5*32.2 feet/sec”2 
All other initial conditions are zero. Figure 15 is a plot of 
the minimum range between the missile and target. The minimum 
distance was 1.13 feet. This was as expected due to the 
changing target trajectory. Figure 16 is the plot of the 
missile/target trajectory in the X-Y plane and Figure 17 is 
the missile/target trajectory in the X-Z plane. Both plots 
demonstrate the ability of the missile using the proportional 
navigation guidance law, to effectively intercept the target. 
Figure 18 is a plot of the missile acceleration vs time. As 
expected, in contrast to the constant velocity simulation, the 
initial portion of the engagement is characterized by a low 
acceleration magnitude. As the engagement progress, however, 
due to the missile having to pull more g's in order to 
compensate for the more rapid trajectory of the target, the 


missile acceleration increased rapidly. 
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Figure 19 is a plot of the targets acceleration vs time. 
Figures 20 and 21 are the seeker head angle rates in the yaw 
and pitch planes, respectively. As with the acceleration of 
the missile, similar characteristics were realized for these 
parameters due to the fact that the missile is constantly 
changing its flight path to compensate for the acceleration of 


the target. 





Figure 15. Range vs Time 
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Figure 18. Missile Acceleration vs Time 





Figure 19. Target Acceleration vs Time 
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MISSILE-TARGET TRAJECTORY IN THE X-Y PLANE 





Figure 16. Trajectory in the X-Y Plane 


MISSILE-TARGET TRAJECTORY IN THB X-Z PLANE 





Figure 17. Trajectory in the X-Z Plane 
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Figure 20. Seeker Head Angle Rate in the Yaw Plane 


REIN Na PCIE JAI LA va FIE 





Figure 21. Seeker Head Angle Rate in the Pitch Plane 
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IV. LUENBERGER OBSERVER 


A. INTRODUCTION 

The purpose of a Luenberger observer for the proportional 
navigation model developed in Chapter III, is to estimate the 
variables in the missile/target engagement that might not be 
measured directly or to reduce noise on the states that are 
measured by the missile. It is not practical to design a 
missile to determine all inputs. Such inputs could be evasive 
maneuvers of the target, initial conditions, noise level, or 
the jamming techniques employed. The use of a Luenberger 


observer is beneficial in these situations. 


B. OBSERVER DESIGN 

The seeker head angle rate is assumed to be a good 
approximation of the line-of-sight rate. The observer will be 
designed to estimate the states of the seeker head, since the 
actual line-of-sight rate of the seeker head can be 
contaminated by noise or jamming. The following equations are 


described in Reference 4. Given a system 
x(k + 1) = Óx(k) + Au(k) (4.1) 
and the output equation 
y(k) =Cx(k) + Du(k) (4.2) 


The observer will approximate x(k) with inputs y(k) and u(k). 
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The observer equation is given as 


2L(k + 1) = FX(k) + Gy(k) + Hulk) (4.3) 


with xhat(k) being the estimate of the original system. The 
error between the state x(k) and the estimated state can be 


defined as 


e(k + 1) = x(k + 1) - X(k + 1) (4.4) 


e(k +1) = ®x(k) + Au(k) - FR(k) - Gy(k) - Hulk) (4.5) 
which results in the error equation of the form 
e(k +1) = (®-GC)x(k) - FR(k) + (A-GD-H)u(k) (4-6) 
If we select F as 
F = ©-GC (4.7) 
Equation (4.4) becomes a state equation of the following form 
e(k + 1) = (©-GO)e(k) + (A-GD-H)u(k) (4.8) 
By selecting H as 
H = A-GD (4.9) 


Equation (4.6), with input u(k), results in a state equation 


without input, of the following form 


e(k +1) = (®-GC) e(k) (4.10) 


If the original system is observable, by selecting the 
observer gain G, we can design the observer such that the 


error approaches zero at a rate that is desired based on the 


design characteristics. 
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The unforced system 


x(k + 1) =0 x (4.11) 


with the observation equation 
y(k) = C x(k) (41-20 


is observable if and only if the rank (N) of the observability 


test matrix 


N= (C "OO Gre (GB) 1C") (4.13) 


is equal to k, the order of the system [Ref.5]. 

The observer design begins by converting the state 
equation for the seeker head, Equation (3.27), from continuous 
time to discrete time representation. With a sampling time of 


0.05 and an unknown disturbance w(k) the equation is given as 


B(k+1) a canus | Bk) ss 
B(k+1) -3.0327 0.3033] |B (k) 3.0327 





0.0013 
alk) z | w(k) (4.14) 


with the output equation measuring the seeker head angle rate 


given as 


M 
k) = [0 (4.15) 
y(k) =f TEN 


The system is observable, so the observer gain can be 
determined. It is usually desirable to select the eigenvalues 
of the observer so that the error between the state and 
estimated state will rapidly become small and will not be 
highly responsive to noise. Because the system is observable, 


the eigenvalues of the observer can be arbitrarily positioned. 
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If the observer gain is high, the observer will quickly drive 
the initial condition error to zero but will make the observer 
more susceptible to noise. By choosing a low observer gain, 
the observer response will be slower but the observer will be 
less affected by noise or jamming. Previous studies have shown 
that a good compromise between the two would be to choose 
continuous time poles of p, and p, to be between the values of 
8.0 to 12.0. Continuous time poles of 10.5 were chosen for the 


observer which results in the discrete time poles given as 


Pewee = oe ee OOO LS (4.16) 


The observer gain matrix is then given as 


0.0299 





The F matrix is equal to (9 - GC) and is given as 


0.9098 a 


-3 10932 1003 83 3 


0.0299 








The H matrix is equal to (A - GD) and with D equal to zero the 
H matrix reduces to 4. The observer equation for the seeker 


head is given as 


^ 0.90 .0297LA -0.00 0.09 
Bien]. 9098 0.029 Bco) >| 0 a yio -| 0902 


4.19 
-3.0327 0.3094 0.0299 3.0327; 7€ ( ) 
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C. SIMULATION RESULTS 
The two previous scenarios given in Chapter III were used 
to test the observer design. Random noise was injected into 
the seeker head to determine if the Luenberger Observer would 
reduce the noise of the states that were directly measured by 
the seeker head and generate the required signals necessary 
for the autopilot. 
1. Constant Velocity Target 
The same initial conditions presented in Chapter III 
for the constant velocity engagement will be used for 
comparison in this section. Figures 22 through 27 are plots of 
the constant velocity engagement without the Luenberger 
Observer but with a random disturbance. Figures 22 and 23 are 
the plots of the trajectories in the X-Y and X-Z planes, 
respectively. The miss distance was 15.6 feet. The trajectory 
in the Y-Z plane depicts how the noise input to the seeker 
head adversely effects the missiles trajectory. Figures 24 and 
25 are the seeker head angle rates in the yaw and pitch plane, 
respectively and Figure 26 is the plot of the commanded 
acceleration as a result of the seeker head angle rates being 
subjected to noise. Figure 27 is the plot of the random noise 


generated by the program. 
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MISSILE-TARGET TRAJECTORY IN THE X-Y PLANB 


MISSILE-TARGET TRAJECTORY IN THE Y-Z PLANE 





Figure 23. Trajectory in the Y-Z Plane with Noise 
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Figure 24. Seeker Head Angle Rate in the Yaw Plane 
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Figure 25. Seeker Head Angle Rate in the Pitch Plane 
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Figure 27. Random Noise Generated vs Time 
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The following plots depict the results with the 
observer included in the simulation with the same initial 
conditions. Figures 28 and 29 are the trajectories in the X-Y 
and Y-Z planes, respectively. In contrast with the Y-Z plane 
plot without the observer, the missiles trajectory is less 
effected by the random disturbance. A miss distance of 1.56 
feet occurred. Figures 30 and 31 are the plots of the seeker 
head angle rates in the yaw and pitch plane, respectively. 
When compared to the angle rate plots of Chapter III, it is 
clear that the observer has reduced the effect of noise to the 
point that the plots are very similar. Figures 32 and 33 are 
the plots of the commanded acceleration and the random noise 
generated, respectfully. The commanded acceleration is also 
less effected by the random disturbance when the observer is 


made part of the missile. 
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MISSILE-TARGET TRAJECTORY IN THE X-Y PLANE 





Figure 28. Trajectory in the X-Y Plane with Observer 


MISSILE- TARGET TRAJECTORY IN THE Y-Z PLANE 





Figure 29. Trajectory in the Y-Z Plane with Observer 
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Figure 30. 


Figure 31. 





Seeker Head Angle Rate in the Yaw Plane 
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Seeker Head Angle Rate 


42 





in the Pitch Plane 








Figure 33. Random Noise Generated vs Time 
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2. Constant Acceleration Target 

The results for the constant acceleration engagement 
were as expected, with the observer reducing the effects of 
noise so that intercept occurred. Figures 34 through 39 are 
plots of the constant acceleration model without the observer. 
Figures 34 and 35 are plots of the trajectories in the X-Y and 
X-Z planes, respectively. These plots do not show the 
variations in the trajectories due to the magnitudes of the 
scales. Figures 36 and 37 are plots of the seeker head angle 
rates in the yaw and pitch planes, respectively. When compared 
to the plots in Chapter III, it is obvious that the missiles 
performance is degraded. A miss distance of 26.5 feet 
occurred. Figures 38 and 39 are plots of the commanded 


acceleration and the random noise generated. 
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MISSILE-TARGET TRAJECTORY IN THE X-Y PLANE 
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Figure 34. Trajectory in the X-Y Plane with Noise 


MISSILE- TARGET TRAJECTOR Y IN THE X-Z PLANE 





Figure 35. Trajectory in the X-Z Plane with Noise 


45 


MIE MA) NAY ANA KWA vs TOR 





Figure 36. Seeker Head Angle Rate in the Yaw Plane 





Figure 37. Seeker Head Angle Rate in the Pitch Plane 
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Figure 39. Random Noise Generated vs Time 
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The following plots depict the constant acceleration 
engagement with the Luenberger Observer. The observer reduced 
the noise so that intercept occurred with a miss distance of 
4.5 feet. The plots in the X-Y and X-Z planes are not included 
due to the magnitude of the scales. Figures 40 and 41 are the 
plots of the seeker head angle rates in the yaw and pitch 
planes, respectfully. When compared with the constant 
acceleration model in Chapter III it is obvious that the 
observer is beneficial in reducing the random disturbance 
added to the seeker head. Figure 42 is the commanded 
acceleration of the missile. 

In the present day environment, with fast maneuvering 
targets, the simulations in this chapter demonstrate the need 
for the observer concept to be included in the missile design 


and also in computer simulation studies. 
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Figure 40. Seeker Head Angle Rate in the Yaw Plane 
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Figure 41. Seeker head Angle Rate in the Pitch Plane 
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Figure 42. Commanded Acceleration with Observer 


The computer program used to generate the plots in this 
chapter is included in the appendix, of which a portion is the 
work completed previously by LT F. Lukenbill in our missile 


and avionics class. 
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V. GROUND TARGET OBSERVER 


A. INTRODUCTION 

Air defense is becoming significantly more complex due to 
the increased speed and manueverability of hostile aircraft 
and the enemy's use of electronic countermeasures (ECM). 
Attacking aircraft will employ jamming as well as other forms 
of ECM to degrade the tracking performance of surface-to-air 
missiles. Therefore, many modern surface-to-air missile 
systems have target observers that track incoming aircraft 
from the ground and uplink the deviations in target position 
and velocity to the missile. 

In this chapter, a simplified Kalman Tracking Filter is 
developed. The tracker consists of velocity and position 
estimation. Once the missile is launched, target deviations 
are computed by subtracting the target observer estimates from 
the current radar estimates. At specified intervals in time, 
the errors in position and velocity are uplinked to the 
missile. After reception, the missile converts the errors into 
the missile body reference frame via the attitude Euler angles 
and then into the line-of-sight frame using the seeker head 
angle transformation. The estimated angular rate deviations 
are then added to the missile’s current estimate and the 


appropriate acceleration commands are generated. 
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B. KALMAN FILTER 

The purpose of the Kalman Filter is to keep track of the 
states of a system given a sequence of noisy measurements. 
This is accomplished by recursively updating an estimate of 
the state by processing a sequence of noisy observations so 
that the effect of the measurement errors are reduced as much 
as possible. 

First, previous estimates of the state x, and the 
covariance matrix P, are extrapolated one step ahead based on 
the assumed systems dynamics. These values are used to compute 
a set of optimum weights or Kalman Gains. The gains are 
applied to the prediction and to a new observation which 
provides an updated estimate of the state and the covariance 
matrix. 

It is assumed that the target tracker has no a priori 
information on the target's position or velocity. Also, the 
error covariance matrix will be set sufficiently large in 
order that the filter gain will be able to make the estimate 
converge to the true value and make the simulation independent 
of any initial estimate states. The random forcing function W, 
is zero mean with covariance Q, and the measurement noise V, is 
zero mean with covariance R. The Kalman Filter algorithm is 
given in Equations (5.1) through (5.5) in the order that the 
filter is implemented in the computer program. The algorithm 


is given as 
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Gy ki P(klk-1) HT (H Pelos H? + R) = 


(5:1) 

RK pe = A (kik-1) CZ = 22 Z (klk-1)) (5.2) 
Pijo 5 (T- GH) Ps) (5.3) 
Rii = 9 Rix (5.4) 

Pikalo = © Pup, 9 * 0 (5.5) 


From Reference 2, the following covariance matrices are 
defined. The purpose of the random forcing function covariance 
matrix Q, is to account for model inaccuracies or for a target 


that has maneuvered. The covariance matrix Q is given as 


0% 0 O| heo o oœ 
Q=|0 o O|=| 0 160? O (5.6) 

0 0 a; 0 0 160? 
Since the actual range of the target is unknown, an estimate 
of the range is used in determining R. The covariance of 
measurement error matrix will then be based on a three degree 
field of view radar, observing the target approaching an area 


20,000 feet from the ground based tracker. 
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$500? 0 0 
Ri == XN 0 (507) 
0 IAS 


The measurement equation y, is given as 


y(k) = H x(k) + v(k) (5.8) 


where H is the measurement matrix, a linear function of the 


state vector, and v(k) the measurement noise. 


C. UPLINK 

The uplink from the ground target observer consists of 
deviations in target position and velocity. The deviations are 
calculated by determining the difference between the actual 


noisy measurement and the predicted state estimate. 


y(k) = H £k (5.9) 


Upon reception, the missile transforms the deviations into the 


missile body frame. 


Vax cos (8)cos(y) cos(8)sin(p) -sin(8)] |Vx 
Vzy| = -sin(w) cos (y) 0 Vy| (5.10) 
Vaz] [Sin(8)cos(w) sin(8)sin(y) cos(6) | |v, 


where Vy, Vy, and V; are the velocity deviations in the x, y, 
and z directions, respecttively, and Vy, Vg,, and Vg, are the 
transformed velocities in the missile reference frame. The 


range deviations are calculated the same and are defined* 
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Fm| [cos(8)cos(y) cos(8)sin(y) -sin(8)]|Fx 
Rgy| = -sin(y) cos (y) 0 Ry| (5.11) 
Raz sin(0)cos(y) sin(0)sin(y) cos(0) Re 


where Ry, R,, and R; are the range deviations and Ry, R,, and 
Rgz are the transformed range deviations in the x, y, and z 
directions, respectively. The transformation to the line-of- 
sight reference frame is accomplished by using the seeker head 


angle transformation and is given as 


Vex} |EOS(0Op)cos(oay) cos(a,)sin(a,) -sin(o,)||V,, 
Val a -sin(0,) cos (a y) 0 Va (5.12) 
Vaz sin(o,)cos(o,) sin(o,)sin(o,) cos(o,) ||Vg 


where c, and c, are the line-of-sight angles in the pitch and 
yaw plane, respectively. Vs, Vsy, and V,, are the transformed 
velocities. The range is transformed using the same 
transformation matrix and is given as 

Rsz| |cos(oe,)cos(o,) cos(o,) sin(a,) -sin(oj)| |Rg, 

Re a -sin(a,) cos (ay) 0 Rey (5.13) 

Ros sin(o,)cos(o,) sin(o,)sin(o,) cos(o,) | |Rs, 
where Ry, Ry, and R, are the transformed ranges in the x, y, 
and z directions, respectively. The final process in the 
uplink transformation is to generate the angular rate 
deviations. These are calculated by dividing the velocity 
deviations by the range and the magnitude of the range 
deviations. The line-of-sight rate in the yaw plane is given 


as Equation (5.14). 
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yaw 7 R+ R, (5215) 


where Ry is given as 


Rp = VU Rox + Rey + Rsz) (5.15) 
and the line-of-sight rate in the pitch plane defined as 


RE 
pitch R+ R, (5.16) 


D. SIMULATION OF THE TARGET OBSERVER 

The constant velocity target in a noisy environment will 
be the first test case for the ground target observer. Figures 
43 through 51 summarize the extensive estimator performance 
calculations for the constant velocity scenario. Figures 43 
through 48 are the histories of the six elements of the state 
vector. An inspection of these figures indicates that the 
estimator tracks the system in position and velocity quite 
well. Figures 49 through 51 are examplese of the Kalman Filter 
Gains vs time. Figures 52 through 54 are the response of the 
missile to the uplink commands. The deviations in position and 
velocity were uplinked at one second intervals. The miss 


distance parameter was observed to be 0.44 feet. 
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TARGET POSITION m the X DIRECTION ws TIME 


** ACTUAL TARGET STATE 


++ KALMAN FILTER ESTIMATE 





Figure 43. Target Position in the X Direction 


TARGET VELOCITY in the X DIRECTION w TIME 


°° ACTUAL TARGET STATS 


++ KALMAN FILTER ESTIMATE . 





Figure 44. Target Velocity in the X Direction 
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TARGET POSITION m the Y DIRECTION vs TIME 


°° ACTUAL TARGET STATE 


++ KALMAN FILTER ESTIMATE . 





Figure 45. Target Position in the Y Direction 


TARGET VELOCITY ia the Y DIRECTION « TIME 


** ACTUAL TARGET STATE 


++ KALMAN FILTER ESTIMATE 





Figure 46. Target Velocity in the Y Direction 
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TARGET POSITION a ue Z DIRECTION s TIME 


** ACTUAL TARGET STATE 


++ KALMAN FILTER ESTIMATE 


TARGET VELOCITY ia the Z DIRECTION w TIME 


** ACTUAL TARGET STATE 


++ KALMAN FILTER ESTIMATE 





Figure 48. Target Velocity in the Z Direction 
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KALMAN GAINS u che X DIRECTION 


*** (POSITION GAIN) 


+++ (VELOCITY GAIN) 


KALMAN GAINS in the Y DIRECTION 


*** (POSITION GAIN) 


*** (VELOCITY GAIN) 





Figure 50. Kalman Gains in the Y Direction 
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KALMAN GAINS ia the Z DIRECTION 


*** (POSITION GAIN) 


*** (VELOCITY GAIN) 





Figure 51. Kalman Gains in the 2 Direction 
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Figure 52. Seeker Head Angle Rate in the Yaw Plane 
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SEEKER HEAD PITCH ANGLE RATE w TIME 





Figure 53. Seeker Head Angle Rate in the Pitch Plane 





Figure 54. Acceleration Commanded 
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For the case Of a maneuvering target, the constant 
acceleration model was used to determine filter performance. 
As expected, filter performance was marginal. Filter 
convergence to an accuracy attained in the constant velocity 
target did not occur. Figures 55 through 63 summarize the 
estimator performance for the constant acceleration scenario. 
Figures 55 through 60 are the six elements of the state 
vector. As shown in these plots, the estimator does not track 
the maneuvering target as well as the constant velocity model. 
Figures 61 through 63 are the Kalman Filter Gains generated 
and Figures 64 through 66 are the missile's response to the 
uplink commands. The uplinks were transmitted using the same 
time interval as in the constant velocity scenario. A miss 


distance of 8.6 feet occurred. 
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TARGET POSITION wm the X DIRECTION w TIME 


** ACTUAL TARGET STATE 


* * KALMAN FILTER ESTIMATE 





Figure 55. Target Position in the X Direction 


TARGET VELOCITY im the X DIRECTION w TIME 


** ACTUAL TARGET STATE 


tt KALMAN FILTER ESTIMATE 





Figure 56. Target Velocity in the X Direction 
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TARGET POSITION m the Y DIRECTION w TIME 


°° ACTUAL TARGET STATE 


+t KALMAN FILTER ESTIMATE 
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* * KALMAN FILTER ESTIMATE 





Figure 58. Target Velocity in the Y Direction 
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TARGET POSITION wa the Z DIRECTION w TIME 


** ACTUAL TARGET STATE 


++ KALMAN FILTER ESTIMATE 


TARGET VELOCITY in the Z DIRECTION w TIME 


** ACTUAL TARGET STATB 


++ KALMAN FILTER ESTIMATE 





Figure 60. Target Velocity in the Z Direction 
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KALMAN GAINS a the X DIRECTION 


eee (POSITION GAIN) 


+++ (VELOCITY GAIN) 


KALMAN GAINS in the Y DIRECTION 


*** (POSITION GAIN) 


*** (VELOCITY GAIN) 





Figure 62. Kalman Gains in the Y Direction 
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KALMAN GAINS a the Z DIRECTION 


*** (POSITION GAIN) 


+++ (VELOCITY GAIN) 





Figure 64. Seeker Head Angle Rate in the Yaw Plane 
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SEEKER HEAD PITCH ANGLE RATE w TIME 


ACCELERATION COMMANDED 





Figure 66. Commanded Acceleration 
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The computer program used to generate the plots in this 


chapter is included in the appendix. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

The use of a three-dimensional model in simulation studies 
is beneficial in that the trajectories of the missile and 
target can be displayed in all three axes, therefore, a better 
understanding of the behavior of the systems involved. 

The simulation studies in Chapter III demonstrated that 
the Proportional Navigation guidance law effectively achieves 
intercept of the target for both the maneuvering and non- 
manuvering target. Chapter IV demonstrated the need for the 
observer concept to be utilized in missile design to reduce 
the effects of noise. Chapter V showed an additional means of 
countering jamming or ECM by employing the ground target 
observer concept. This demonstrated that the uplinks in target 
position and velocity, for the constant velocity model, were 


beneficial in a noisy environment. 


B. RECOMMENDATIONS 

Further studies could include a Kalman Filter in the 
missile design instead of the Luenberger Observer. It is 
recommended that future studies include the use of a two part 
Kalman Filter with the ground target observer. One that 


estimates the states and the other being a maneuver detector. 
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The results could be compared with the present filter along 


with ECM effects to determine the effectiveness of the two. 
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APPENDIX 


A. MISSILE-TARGET ENGAGEMENT PROGRAM 
% 
% 
% 
% an asymptotic observer. 
% DEFINE STATES 
% 
5 Missile 
% 
% ms=[xm - missile position 
5 xdm - missile velocity 
$ ym - missile position 
5 ydm - missile velocity 
% zm - missile position 
% zdm - missile velocity 
% 
Am=[0 10000 
000000 
000100 
000000 
00000 1 
o 0o 0 0 0 0j]; 
Bm=[0 O O 
100 
000 
010 
0 0 O 
O O 1]; 
% 
% Seeker Head 
% 
% betap=(beta pitch 
% 
% betay=(beta yaw 
% 
% 
% 


Asp-[O 1;-100 -20.0]; 
Asy-[O 1;-100 -20.0]; 
Bsp=[0;100)]; 


in 
in 
in 
in 
in 
in 


N N NG NG XX 


This program simulates a 3 dimensional target/missile 
engagement using basic proportional navigation with 


direction 
direction 
direction 
direction 
direction 
direction] 


- seeker head pitch angle 


betad pitch - seeker head pitch angle rate 


- seeker head yaw angle 
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betad yaw - seeker head yaw angle rate] 


sy=[0;100]; 
Autopilot 


gammamp=(gammam pitch - pitch body angle 

gammadm pitch - pitch body angle rate |) 
gammamy=(gammam yaw - yaw body angle 

gammadm yaw - yaw body angle rate] 


de AP AL ge oe ge oe de Ud 


Ap=(-3 0;0 -3]; 


Bp=[3 0;0 3]; 


% 
% Target 
% 
% ts=(xt - target position in x direction 
$ xdt - target velocity in x direction 
% yt - target position in y direction 
% ydt - target velocity in y direction 
% ZE - target position in z direction 
% zdt - target velocity in z direction] 
% 
At=(0 I0 000 

000000 

000100 

000000 

00000 1 

oo a favor, 
Bt=(0 0 O 

100 

000 

010 

000 

O O 1); 
$ 
% DISCRETE REPRESENTATION 
% 
dat=.05; 


(phisp,delsp]= c2d(Asp,Bsp,dt) ; 
(phisy,delsy]= c2d(Asy,Bsy,dt) ; 
(phim,delm]= c2d(Am,Bm,dt); 
(phia,dela]= c2d(Ap,Bp,dt) ; 
(phit,delt]= c2d(At,Bt,dt); 


tfinal= 20.0; 
kmax=tfinal/dt +1; 


% INITIAL VALUES 
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+ 
t navigation ratio 
% 

NR=(4.0 0;0 4.0]; 


% initial seeker head angles and angle rates 
$ 

beta_pitch0 =0.0; 

betad_pitch0=0.0; 

beta yawo =0.0; 

betad yawO =0.0; 


betap(:,l)=[beta pitcho 
betad pitchO ]; 

betay(:,1)-[beta yawO 
betad yawO J]; 


$ initial observer angles and angle rates 


betae pitch0-0.0; 

betade_pitch0=0.0; 

betae yaw0-0.0; 

betaed yaw0=0.0; 

$ 

t 

betaep(:,1)-[betae pitchO 
betade pitch0], 

betaey(:,1)-[betae yawO 
betaed yaw0]; 

4 

x 

% initial missile body angle rates 

$ 

gammadm pitcho=o.o; 

gammadm yaw0 =o.o; 


gammadm(:,1)=(gammadm_pitcho 
gammadm yawO ]j; 


$ 

% initial states for missile flight profile 

$ 

xm0 =000.0; 

xdm0=2000.00; 

ymO =000.0; 

ydm0=00; 

zmO -00; 

zdm0=00; 

ms (:,1)=(xmo0 
xdmo 
ymo 


75 


ydmo 


zmo 
zdm0]; 
3 
$ initial states for target flight profile 
% 
XTO -20000; 
xdt020000; 


yt0 =050.0; 
ydt0=0300.0; 
ztO =0051.0; 
zdt0=000.0; 


ts(:,1)=(xto0 
xdto 
yto 
yato 
zto 
zdt0o); 


t 

% initial range information 
% 

rxO-(xtO-xm0); 

rx(1)=rx0; 

ry0=(yt0-ymo0) ; 

ry (1) =ry0; 

rzO=(ztO-zm0); 

rz(l)=rzo; 

rmO-sqrt(xm0O^2 + ymO“2 + Zm0^2); 
rm(1)-rm0; 

rtO-sqrt(xt0^2 * ytO*2 + zt0%*2); 
rt(1)=rto; 

ro =sqrt((xtO-xm0)*2 + (ytO-ym0)*2 + (ztO-zm0) 2); 
r(1)=r0; 


d15tz[0.0013:0.0513]: 
rand('normal'); 


2 

% initial time 
$ 

time(1)=0.0; 

$ 

X SIMULATION 

$ 


for (i=1:kmax-1) 
% missile and target flight path angles 
gammam _ pitch(i)=atan2(ms(6,i) ,sqrt(ms(2,i)*2+ms(4,i)*2)); 


gammam vaw(i)-atan2(ms(4,i),ms(2,1)); 
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gammat pitch(i)-atan2(ts(6,i),sqrt(ts(2,i)^2*ts(4,1)^2)); 
gammat yaw(i)=atan2(ts(4,1i),ts(2,1)); 


i target velocity (magnitude) 


vt (i)=sqrt(ts(2,i)*2+ts(4,i)*2+ts(6,i)*2); 
vm (i)=ssqrt(ms(2,1)*2+ms(4,1)*2+ms(6,1)*2); 
& 
t line-of-sight angles 
% 


sigma pitch(i)-atan2((ts(5,i)-ms(5,i)),sqrt((ts(1,i)-ms(1,i) 
“2... 


) 
+(ts(3,1)-ms(3,1))*2)); 


sigma yaw(i)-atan2((ts(3,i)-ms(3,i)),((ts(1,1i)-ms(1,1)))); 


sigma(:,i)=(sigma pitch(i) 
sigma yaw(i) ]; 


% update seeker head 
$ 
wk(1)=rand*.9; 


betap ( : ,1+1)=phisp*t*betap(: , i) + 
delsp*sigma_pitch(1)+dist*wk(i); 


betay(:,i*l)-phisy*betay(:,i) + delsy*sigma yaw(i)+dist*wk(i); 


$ 
% set up observer estimates 
% 

y0(1)=0.0; 

y1(1)=0.0; 

% 


% 

gain=(-0.0031;0.0299]; 

c=(0 1]; 

f=phisp-gain*c; 

y0(i+1)=c*betap(:,i); 

betaep(:,i+1)=f*betaep(:,1i) + gain*yO(i) ~ 
delsp*sigma pitch(i); 


y1(i+1)=c*betay(:,i); 

fl-phisy-gain*c; 

betaey(:,i*1)-fl*betaey(:,i) + gain*y1(i) + 
delsy*sigma yaw(i); 


n 


IP IP AP AP 


9e 9e 


9e 9e 9e 9e 


seeker head angle rate vector 
betad(:,i)s(betaep(2,i) 
betaey(2,i)]: 

update autopilot 
acceleration commanded = vm(i) * gammadm(i) 
gammadm(i) = Nav constant * line of sight rate 
gammadm(:,i+1)=phia*gammadm(:,i) + dela*NR*betad(:,i); 
vm_pitch(i)=vm(i)*cos(gammam _yaw(i)-sigma_yaw(i)); 
vm yaw(i)-vm(i)*cos(gammam pitch(i)); 
limit commanded acceleration to approximately 20 g's 

if (vm_yaw(i)*gammadm(2,i))<=640.0 

amc y(i)-(vm yaw(i)*gammadm(2,i)); 
else 


amc y(i)-2640.0; 
end 


if(vm pitch(i)*gammadm(1,i))«2640.0 
amc p(i)-2(vm pitch(i)*gammadm(1,1i)); 
else 
amc p(i)2640.0; 
end 
magnitude of commanded acceleration 
acom(i)-sqrt(amc y(i)^2 -* amc p(i)^2); 
missile acceleration vector components 
xddm(i)-(gammadm(1,i)*ms(6,i)*cos(gammam yaw(i))); 
yddm(i)-amc y(i)*cos(gammam yaw(i)); 
zddm(i)-amc p(i)*cos(gammam pitch(i)); 
total missile acceleration magnitude 


am(i)=sqrt(xddm(i)*2 + zddm(i)*2 + yddm(i)*2 ); 


missile input vector 
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de de ap 9e AP ge 


de ap ap de oe oe 


de AP ap oe de oe oe 


de Je X oe 


um=(xddm (1) 
yddm (1) 
zddm(i)); 

update missile 


ms (:,1+1)=phim*ms(:,i) + delm*un; 


set target acceleration components 


if (r(i)<=21000.0) 


xddt (i)=-5.0*32.2*sin(sigma_yaw(i)); 
yddt(i)s» 3.5*32.2*cos(sigma yaw(i)); 
zddt(i)= 3.5*32.2*cos(sigma_pitch(i)); 


else 


xddt(1)=0.0; 
yddt(1)=0.0; 
zddt(i)=0.0; 


end 


target acceleration magnitude 
at(i)=sqrt(xddt(i)*2 + yddt(i)*2 + zddt(i)*2); 
target input vector 
ut=[xddt(i) 
yddt(i) 
zddt(i)]; 
update target 


ts(:,i+1)=phit*ts(:,i) + delt*ut; 


missile and target trajectory data 


missile(i,:)=(ms(1,i) ms(3,i) ms(5,i)); 
target(i,:) =(ts(1,i) ts(3,i) ts(5,1)); 


update range information 


r(i+l)=sqrt((ts(1,1+1)-ms(1,1+1))*2 + 


(ts(3,1+1)-ms(3,1+1))*2... 


| + (ts(5,i+1)-ms(5,i+1))^2); 
rm(1+1)=sqrt(ms(1,i+1)^2 +ms(3,i+1)^2 +ms(5,i+1)^2); 
rt(i+1)=sqrt(ts(1,i+1)^2 +ts(3,i+1)^2 +ts(5,i+1)^2); 
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$ update time vector 
$ 
time(i+1)=time(i) + dt; 
% 
% check to see if engagement is at closest point of 
approach 
if(r(i)<r(i+1)) break,end 
end; 
% 
X PLOT RESULTS 
% 
% range data 
% 


plot(time,r),grid,xlabel(”TIME/),ylabel(*FEET/) 
title(’RANGE VS TIME’) 
text( 7,20000,['rng-s',num2str(r(i)),' feet']); 
meta msltgt,pause,clg 


subplot(211),plot(time,rm),grid,xlabel('TIME'/),ylabel('FEET') 
title('MISSILE RANGE VS TIME'); 
subplot(212),plot(time,rt),grid,xlabel('TIME'/),ylabel('FEET') 
title(’TARGET RANGE VS TIME’); 

meta jerryl,pause,clg 


% 

plot(ms(1,:),ms(3,:),ts(1,:),ES(3,99 grad 
title('MISSILE-TARGET TRAJECTORY IN THE X-Y PLANE") 
xlabel('X AXIS IN FEET'),ylabel('Y AXIS IN FEET’) 
gtext('MISSILE') 

gtext(’TARGET’ ) 

meta msltgt,pause,clg 


plot(ms(1,:);ms(5,2), tsi.) ES(S ON GRE 
title('MISSILE-TARGET TRAJECTORY IN THE X-Z PLANE’) 
xlabel(’X AXIS IN FEET’) ,ylabel(’Z AXIS IN FEET’) 
gtext('MISSILE') 

gtext ('TARGET') 

meta msltgt,pause,clg 


% missile and target velocity 

% 

plot(time(1:i) ,vm) ,grid,xlabel('TIME'),ylabel('FEET/SEC') 
title(’MISSILE VELOCITY VS TIME’); 

meta jerryl,pause,clg 


plot(time(1:i),vt),grid,xlabel(’TIME’) ,ylabel(’FEET/SEC’ ) 
title(’TARGET VELOCITY VS TIME’); 
meta jerryl,pause,clg 


% missile and target acceleration 
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plot(time(1:i),am),grid,xlabel('time'),ylabel('feet/sec^2') 
title(’Missile acceleration vs Time’) 

meta jerryl,pause,clg 

plot(time(1:i) ,at) ,grid,xlabel('time') ,ylabel('feet/sec*2') 
title(’Target acceleration vs Time’) 

meta jerryl,pause,clg 


% seeker head angle rates 


plot(time(1:i-1),betay(2,(1:i-1))),grid,xlabel('time'),ylabe 
l('radians') 

title(’Seeker head Yaw angle rate vs Time’) 

meta jerryl,pause,clg 

plot(time(1:i-1) ,betap(2,(1:1-1))),grid,xlabel('time') ,ylabe 
l('radians') 

title(’Seeker head Pitch angle rate vs Time’) 

meta jerryl,pause,clg 


% commanded acceleration 


plot(time(1:i),acom) ,grid,title(’acceleration commanded’ ) 
meta jerryl,pause,clg 


END 


81 


B. MISSILE-TARGET PROGRAM WITH TARGET OBSERVER 


clear 

erase msltgt.met 
lerase jerryl.met 
clg 


position 
velocity 
position 
velocity 
position 
velocity 


in 
in 
in 
in 
in 
in 


N NG NG XX 


This program simulates a 3 dimensional target/missile 
engagement using basic proportional navigation with 
an Luenberger observer and a ground target observer. 


direction 
direction 
direction 
direction 
direction 
direction] 


- seeker head pitch angle 


pitch - seeker head pitch angle rate 


- Seeker head yaw angle 


yaw - seeker head yaw angle rate] 


$% 
$% 
t 
% 
% 
% 
% DEFINE STATES 
t 
% Missile 
% 
% ms=[xm - missile 
4 xdm - missile 
% ym - missile 
% ydm - missile 
$ zm - missile 
% zdm - missile 
% 
Am=(0 10000 
000000 
000100 
000000 
000001 
000000]; 
Bm=(0 0 0 
100 
000 
0106 
000 
O O 1j]; 
% 
t Seeker Head 
% 
% betap=[beta pitch 
t betad 
t betay=(beta yaw 
t betad 
% 
t 


Asp-[O 1;-100 -20.0]; 
Asy=[0 1;-100 -20.0]; 


Bsp=[0;100]; 
Bsy=[0;100]; 
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Je oe AP AP AP AP AP AP 


Ap=(-3 


Autopilot 


gammamp=(gammam pitch - pitch body angle 


gammadm pitch - pitch body angle rate 


gammamy=[gammam yaw - yaw body angle 
gammadm yaw - yaw body angle rate] 


010—317 


Bp=[3 0;0 3]; 


3» doe oe AP AP AP AP AP AP AP 
OoOooooo 


QOOOOO fr 


QDOOOHF O 
OOPOOO 


dt=.05; 


Target 

ts=(xt 
xdt 
yt 
yat 
ZE 
zdt 

0000 

0000 

0100 

0000 

0001 

oOo 0 0 0J; 

0 

0 

0 

0 

0 

1]; 


- target 
- target 
- target 
- target 
- target 
- target 


position 
velocity 
position 
velocity 
position 
velocity 


DISCRETE REPRESENTATION 


(phisp,delsp)= c2d(Asp,Bsp,dt); 
(phisy,delsy)= c2d(Asy,Bsy,dt); 


(phim,delm]- c2d(Am,Bm,dt); 
(phia,dela)= c2d(Ap,Bp,at); 
(phit,delt]s c2d(At,Bt,dt); 


tfinal= 20.0; 
kmax=tfinal/dt +1; 


% 
% 


INITIAL VALUES 


83 


in 
in 
in 
in 
in 
in 


N NKK X X 


direction 
direction 
direction 
direction 
direction 
direction) 


] 


navigation ratio 


% 
$ 
NR=[4.0 0;0 4.0]; 


% initial seeker head angles and angle rates 
% 

beta pitch0 =0. 
betad pitcho=o. 
beta yawo =0. 
betad yawO =0. 


OOOO 


betap(:,1)=[beta pitcho 
betad pitcho ]; 

betay(:,1)-2[beta yawO 
betad yawO |]; 


+ initial line-of-sight angular rates 


sigmad pitch0-0.0, 
sigmad yawO  -0.0; 


sigmad(:,l)=[sigmad pitcho 
sigmad yawo ]; 


x initial observer angles and angle rates 


betae pitch0-0.0, 

betade pitch0=0.0; 

betae yaw0=0.0; 

betaed yawo=o.o; 

$ 

$ 

betaep(:,l)=(betae pitchO 
betade pitch0]; 

betaey(:,l)=[betae yawo 
betaed yaw0]; 


initial missile body angle rates 


IP IP de ge 


gammadm pitch0-0.0; 
gammadm yawO =0.0; 


gammadm(:,1)=(gammadm pitcho 
gammadm yaw0 ], 

% 

% initial missile states 

$ 

xm0 =000.0; 

xdm0=2000.00; 
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ymo =000.0; 


ydm0z00; 

zmO 200; 

zdm0=00; 

ms(:,1)=(xmo0 
xdmO 
ymo 
yamo 
zmo 
zam0]; 

% 

% initial target states 

% 


xtO =20000; 
xdt020000; 
ytO 2050.0; 
ydt0=1000.0; 
ztó =050.0; 
zato=000.0; 


ts(:,l)=[xto 


xdto 

yto 

ydto 

ZtO 

zdt0]; 
4 
4 initial range information 
X 


rxO-(xtO-xm0); 

rx(1)=rx0; 

ryO-(yt0-ym0); 

ry(1)=ry0; 

rz0=(zt0-zm0); 

rz(1)=rz0; 

rmO-sqrt(xm0^2 + ym0”*2 + zm0*2); 
rm(1)=rmo0; 

rto=ssqrt(xt0*2 + yto*2 + zto^2); 
rt(1)=rt0; 

r0 =sqrt ((xt0-xm0)*2 + (yt0-ym0)^2 t (zt0-zm0) 12), 
r(1)=r0; 


% Initial Target Tracker values 

y(:,l)=[0;0;0]; 
xhatl =(0;0;0;0;0;0); 
pkkm1=(10*9 00000 

O 10*9 0000 

00 10*9 00 0 

00 0 10*9 0 0 

0000 109 0 
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00000 10*9]; 


rand('normal'); 


DOO O Oo 


o 
H 
I 
ns 
HOO QO OO 


DO OO FO 
DOH QO 2o 


); 
dist-[(0.0013:0.0513], 


Continuous to Discrete 
(phi ,de1]-c2d(aT,bT,dt), 


initial time 


SIMULATION 
(i=1:kmax-1) 
ES 
% missile and target flight path angles 


% 
% 
% 
% 
% 
time(1)=0.0; 
5 
$ 
J 
f 
J 


gammam pitch(i)-atan2(ms(6,i),sqrt(ms(2,i)^2*ms(4,1)^2)); 
gammam yaw(i)zatan2(ms(4,i),ms(2,i)); 
gammat pitch(i)-atan2(ts(6,i),sqrt(ts(2,i)^2*ts(4,1)^2)); 
gammat_yaw(i)=atan2(ts(4,i),ts(2,i)); 


% target velocity (magnitude) 


vt (i)=ssqrt(ts(2,i)*2+ts(4,i)*2+ts(6,1)*2); 
vm (i)=sqrt(ms(2,i) “2+ms(4,i1) '2+ms(6,1) 2); 


9e AP 


line-of-sight angles 
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% 
sigma pitch(i)-atan2((ts(5,i)-ms(5,i)),sqrt((ts(1,i)-ms(1,i) 
E... 
HOCICO ETE 
sigma yaw(i)-atan2((ts(3,i)-ms(3,i)),((ts(1,i)-ms(1,i)))); 
sigma(:,i)-(sigma pitch(i) 
sigma yaw(i) ]; 
4 update seeker head 


wk (1)=rand*.9; 


betap(:,i+1)=phisp*betap(:,1) + 
delsp*sigma_pitch(i)+dist*wk(i); 

betay (:,1+1)=phisy*betay(: ,i) + 
delsy*sigma_yaw(i)+dist*wk(i); 
% 
% set up observer estimates 
4 


y0(1)=0.0; 
y1(1)=0.0; 
% 


t 

gain=(-0.0031;0.0299]; 

c=(0 1]; 

f=phisp-gain*c; 

y0(1+1)=c*betap(:,i); 

betaep(:,1+1)=f*betaep(:,1) + gain*yO(i) + 
delsp*sigma_pitch(i); 


yl(i*1)sc*betay(:,i); 

fl-phisy-gain*c; 

betaey(:,i+1)=f1*betaey(:,1i) + gain*yl(i) * 
delsy*sigma yaw(i); 


seeker head angle rate vector 
with uplinks at one second intervals 


9e AP AP AP AP 


if (j==40) 
betad(:,i)=[ (betaep(2,i)+sigmad(1,i)) 
(betaey(2,i)+sigmad(2,i))]; 
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de oe 


oe oe oe oe oe 


j=0; 
else 
betad(:,1)=[betaep(2,1) 
betaey(2,1)]; 

end 
update autopilot 
acceleration commanded = vm(i) * gammadm(i) 
gammadm(i) = Nav constant * line of sight rate 
gammadm(:,i*1)-phia*gammadm(:,i) + dela*NR*betad(:,i) ; 
vm _pitch (1) =vm(1) *cos(gammam_ yaw(i)-sigma yaw(i)); 
vm "yaw(i)-vm(i)*cos(gammam pitch(i)); 
limit commanded acceleration to approximately 20 g's 

if (vm yaw(i)*gammadm(2,i))«2640.0 

amc y(i)-2(vm yaw(i)*gammadm(2,1i)); 
else 


amc y(i)2640.0; 
end 


if(vm pitch(i)*gammadm(1,i))«-2640.0 
amc p(i)-(vm pitch(i)*gammadm(1,1i)); 
else 
amc p(i)2640.0; 
end 
magnitude of commanded acceleration 
acom(i)-sqrt(amc y(i)^2 * amc p(i)^2); 
missile acceleration vector components 
xddm(i)z(gammadm(1,i)*ms(6,i)*cos(gammam yaw(i))); 
yddm(i)zamc y(i)*cos(gammam yaw(i)); 
zddm(i)-amc p(i)*cos(gammam pitch(i)); 
total missile acceleration magnitude 


am(i)=sqrt(xddm(i)*2 + zddm(i)*2 + yddm(i)*2 ); 


missile input vector 
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of oe oe 9e AP AP 


AP AP AP ge 9e oe ge 9e JE oe Je JE 


Je Je 9e oe 


9e de oe x 


ums ( xàdm(i) 
yddm(1) 
zddm(i)]; 
update missile 
ms (:,1+1)=phim*ms(:,i) + delm*um; 
set target acceleration components 
if (r(i)<=00000.0) 
xddt(i)= -S5.0*32.2*sin(sigma_yaw(i)); 
yddt(1)= 3.5*32.2*cos(sigma_yaw(i)); 
zddt(i)= 3.5*32.2*cos(sigma pitch(i)); 
else 
xddt(1)=0.0; 
yddt(1)=0.0; 
zddt(1)=0.0; 
end 
target acceleration magnitude 
at(i)=sqrt(xddt(i)*2 + yddt(i)*2 + zddt(i)*2); 
target input vector 
ut=(xddt (i) 
yddt (i) 
zddt(i)]; 
update target 


ts(:,i+1)=phit*ts(:,i) + delt*ut; 


missile and target trajectory data 


missile(i,:)={ms(1,i) ms(3,i) ms(5,i)]; 
target (lan) =a(ts(i1,i)ots(3,1) tg(5,i)); 


Target traker 


Input initial conditions 


es 1 1) /€s(2,1),€9(3,1),€3(4,1);88(5,1);ts(6,1)]; 


Plant noise covariance 


Q1-[160^2 O 0;0 16052 0:0 O 160742], 


39 


9e 9e 9e 9e 9e AP AP de AP 


9e 9e 


Measurement noise covariance 


R=(500*2 0 O 
0 1252 O 
0: 001222]: 

Input equations 
Q=[160*2;160*2;160*2]; 
nusz[sqrt(500)^2;sqrt(125^2);sqrt(12^2)]*rand; 
w(:,1) =sqrt(Q)*rand; 


Update the plant states 


x(:,i*l)sphi*x(:,i) +del*w(:,i); 
y(:,i+1)=h*x(:,i+1) +nu; 


Compute Kalman gains 

G = pkkml*h'/*inv(h*pkkml*h'-R); 
Define gain matrix for plotting 

g(:,1) =[G(1:2,1);613:34,2) CMS SY 
Measurement step 


xhat (:,i+1) =xhat1+(G *(y(:,i+1)-h*xhat1)); 
pkk =(eye(6)-G *h)*pkkm1; 


Movement step 


xhatl=phi*xhat (:,i+l); 
pkkmi=phi*pkk*phi'+ del*Q1*del'; 


compute the range deviation 
xhatr -[xhat(1,i);xhat(3,i);xhat(5,1)]; 
yr  =y(:,i); 
xr(:,i)z (yr - xhatr); 
compute the velocity deviation 
xhatv  -[xhat(2,i);xhat(4,i);xhat(6,1)]; 


yv z[ts(2,1);ts(4,1);ts(6,1)]; 
xv(:,1)= (yv > xhatv); 
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% Transformation to body axes 


all-cos(gammam pitch(i))*cos(gammam yaw(i)); 
al2-cos(gammam pitch(i))*sin(gammam yaw(i)); 
al3=(-sin(gammam pitch(i))); 
a21-(-sin(gammam yaw(i))); 

a22=cos(gamman yaw(i)); 

a23=0.0; 

a3l-sin(gammam pitch(i))*cos(gammam yaw(i)); 
a32-sin(gammam pitch(i))*sin(gammam yaw(i)); 
a33=cos(gammam pitch(i)); 


% velocity transformation 


vbx(i)=xv(1,i)*all+xv(2,i)*a12+xv(3,i)*a13; 
vby(i)=xv(1,i)*a21+xv(2,1)*a22+xv(3,i)*a23; 
vbz(i)=xv(1,1)*a3l+xv(2,1)*a32+xv(3,1)*a33; 


€ range transformation 


rbx(i)=xr(1,i)*all+xr(2,1)*al2+xr(3,1)*al3; 
rby(i)=xr(1,i)*a21+xr(2,i)*a22+xr(3,1)*a23; 
rbz(i)=xr(1,1) *a31+xr(2,1i) *a32+xr(3,1) *a33; 


$ line-of-sight transformation 


sll-cos(sigma pitch(i))*cos(sigma yaw(i)); 
sl2-cos(sigma pitch(i))*sin(sigma yaw(i)); 
sl3-(-sin(sigma pitch(i))); 
s21-(-sin(sigma yaw(i))); 

s22-cos (sigma yaw(i)); 

s23=0.0; 

s3l-sin(sigma pitch(i))*cos(sigma yaw(i 
s32-sin(sigma pitch(i))*sin(sigma yaw(i 
533=c0s (sigma _pitch(i)); 


Naat? eee 


) a; 
) ; 
% velocity transformation 
vsx(i)=vbx(i)*sll+vby(i)*sl2+vbz(i)*s13; 
vsy (1) =vbx(1) *s21+vby (1) *s22+vbz (1) *s23; 
vsz(i) =vbx (1) *s31l+vby (1) *s32+vbz (1) *s33; 
rsx(i)=rbx(i)*sll+rby(i)*sl2+rbz(1)*s13; 
rsy(i)srbx(i)*s21-*rby(i)*s22-*rbz(i)*s23; 
rsz(i)=rbx(1)*s31+rby(1)*s32+rb2(1)*s33; 
% compute the magnitude of the range deviation 


rd(i)=sqrt(rsx(i)*2+rsy(i)*2+rsz(i)%*2); 


% generate the angular rate deviations 


91 


Sigmad pitch(i+l)= vsy(i)/(r(1)+rd(i)); 
sigmad yaw(i*1)- vsz(i)/(r(i)*rd(i)); 


sigmad(:,i+l)=[sigmad pitch(i) 
sigmad yaw(i)]; 


% update range information 
% 


r(i+1)=sqrt((ts(1,i+1)-ms(1,i+1))"2 + 
(ts(3,i+1) -ms(3,1+1))*2... 
+ (ts(5,1+1)-ms(5,1+1))*2); 
rm(i+1)=sqrt(ms(1,1+1)*2 +ms(3,1+1)*2 +ms(5,i+1)"2); 
rt (1+1)=sqrt (ts(1,1+1)*2+ts(3,i1+1)*2+ts(5,i+1)"2); 


4 update time vector 
a 
time(i+1)=time(i) + dt; 
% 
% check to see if engagement is at closest point of 
approach 
% 
if(r(i)<r(i+1)) break,end 
end; 
% 
% PLOT RESULTS 
% 
% range data 
x 


plot(time,r) ,grid,xlabel(’TIME’) ,ylabel(’ FEET’ ) 
title(’RANGE VS TIME’) 
text( 7,20000,['rngs',nun2str(r(i)),' feet']); 
meta msltgt,pause,clg 


subplot(211),plot(time,rm),grid,xlabel('TIME'),ylabel('FEET') 
title(’MISSILE RANGE VS TIME’); 

subplot(212) ,plot(time,rt),grid,xlabel (“TIME”) ,ylabel (“FEET”) 
title('TARGET RANGE VS TIME”); 

meta msltgt,pause,clg 

% 

% 

% missile and target velocity 

% 

$plot(time(1:i),vm),grid,xlabel('/TIME'),ylabel('FEET/SEC') 
$title('MISSILE VELOCITY VS TIME’); 

%meta msltgt,pause,clg 


tplot(time(1:i), vt),grid, xlabel (“TIME”) ,ylabel('FEET/SEC") 


$title(’TARGET VELOCITY VS TIME’); 
tmeta msltgt,pause,clg 
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$% missile and target acceleration 


plot(time(1:i) ,am) ,grid,xlabel('TIME'),ylabel('FEET/SEC*2') 
title('MISSILE ACCELERATION vs TIME’) 

meta msltgt,pause,clg 
plot(time(1:i),at),grid,xlabel('time'),ylabel('feet/sec^2') 
title('TARGET ACCELERATION vs TIME') 

meta msltgt,pause,clg 


% seeker head angle rates 


plot(time(1:i-1),betaey(2,(1:i-1))),grid,xlabel('time'),ylab 
el('radians') 

title(’SEEKER HEAD YAW ANGLE RATE vs TIME’) 

meta msltgt,pause,clg 

plot(time(1:i-1) ,betaep(2,(1:1-1))),grid,xlabel('time') ,ylab 
el('radians") 

title(’SEEKER HEAD PITCH ANGLE RATE vs TIME’) 

meta msltgt,pause,clg 


% commanded acceleration 


plot(time(1:i) ,acom) ,grid,title('ACCELERATION COMMANDED/) 
xlabel('TIME'),ylabel('FEET/SEC^2') 
meta msltgt,pause,clg 


t Kalman estimates 


plot(time(1:i-1),x(1,(1:i-1)),'*',time(1:i-1),xhat(1,(1:i-1) 
MT”) 

xlabel('TIME'),ylabel('FEET'),grid 

title('TARGET POSITION in the X DIRECTION vs TIME’) 
gtext('** ACTUAL TARGET STATE ’) 

gtext(’++ KALMAN FILTER ESTIMATE') 

meta msltgt,pause,clg 


plot(time(1:i-1),x(2,(1:i-1)),'*',time(1:i-1),xhat(2,(1:i-1) 
It") 

xlabel (“TIME”) ,ylabel(“FEET/SEC”) , grid 

title(’TARGET VELOCITY in the X DIRECTION vs TIME’) 
gtext(’** ACTUAL TARGET STATE ^") 

gtext(’++ KALMAN FILTER ESTIMATE’ ) 

meta msltgt,pause,clg 


plot(time(1:1-50),g(1,(1:i-50)),'*',time(1:i-50),g(2,(1:1-50 
+”) 

title('KALMAN GAINS in the X DIRECTION’) ,grid,xlabel(’TIME’), 
gtext('*** (POSITION GAIN) ’) 

gtext(’+++ (VELOCITY GAIN) ’) 

meta msltgt,pause,clg 
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% Kalman estimates 


plot(time(1:1-1),x(3,(X:1-1)),'*',time(l1:1—1), xhae(3 (SEE 
! E) 

«Label ("TIME") ylarel(-FEET"),grid 

title( TARGET POSITION in the Y DIRECTION vs TIME’) 

gtext(’** ACTUAL TARGET STATE  ') 

gtext(’++ KALMAN FILTER ESTIMATE’ ) 

meta msltgt,pause,clg 


plot (time(1:i-1),x(4,(1:i-1)),'*” time(1:i-1),xhat(4,(1:i-1) 
Ja 

xlabel(/TIME'),ylabel('FEET/SEC') grid 

title(/TARGET VELOCITY in the Y DIRECTION vs TIME') 
gtext('** ACTUAL TARGET STATE  ') 

gtext(’++ KALMAN FILTER ESTIMATE’) 

meta msltgt,pause,clg 


plot(time(1:1-50) ,g(3,(1:1-50)),’*’,time(1:1-50) ,g(4, (1:1-50 
ZY 

title(’KALMAN GAINS in the Y DIRECTION’) ,grid,xlabel(’TIME’), 
gtext(’*** (POSITION GAIN) ') 

gtext(’+++ (VELOCITY GAIN) ’) 

meta msltgt,pause,clg 


plot(time(1:i-1),x(5,(1:i-1)),'*',time(1:1i1-1),xhat(5 P] (15m 
yet) 

xlabel('TIME'),ylabel('FEET'),grid 

title(’TARGET POSITION in the Z DIRECTION vs TIME’) 
gtext(’** ACTUAL TARGET STATE  ') 

gtext(’++ KALMAN FILTER ESTIMATE’ ) 

meta msltgt,pause,clg 


plot (time(1:i-1),x(6,(1:i-1)),**” time(1:i-1) xhat(6,(1:1-1) 
jit 

xlabel (“TIME”) ,ylabel(“FEET/SEC”) , grid 

title('TARGET VELOCITY in the Z DIRECTION vs TIME") 
gtext('** ACTUAL TARGET STATE 2 

gtext('++ KALMAN FILTER ESTIMATE’) 

meta msltgt,pause,clg 


plot(time(1:1-50) ,g(5,(1:1-50)),/#/” ,time(1:1-50) ,g(6, (1:1-50 
Jona a 

title(’KALMAN GAINS in the Z DIRECTION’) ,grid,xlabel(’TIME’), 
gtext(’*** (POSITION GAIN) ’) 

gtext(’+++ (VELOCITY GAIN) ’) 

meta msltgtl1,pause,clg 

end 
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